package align_gaps;

require Exporter;


our @ISA = qw(Exporter);
our $VERSION = 1.00;
our @EXPORT = qw(   map_pos_via_aligned 
					map_pos_from_aligned
					map_pos_to_aligned
					expand_prot_with_gapstr
					expand_cdna_with_gapstr
					expand_prot_with_gaps
					expand_cdna_with_gaps
					gapstr_to_gaps
					gaps_to_gapstr
					aligned_seq_to_gaps	);
use strict;
use warnings;
use Carp;


#_________________________________________________________________________________________

#	gapstr_to_gaps

#		input:

#			$gap_str
#			contains gaps in begin_pos:length|begin_pos:length format
#			0:63|119:6|263:1

#		returns:

#			reference to array of [[beg_pos:length][beg_pos:length]]

#	gaps_to_gapstr

#		reverses above

#_________________________________________________________________________________________
sub gapstr_to_gaps($)
{
	return [map {[split '\:', $_]; } split '\|', $_[0]];
}
sub gaps_to_gapstr(\@)
{
	return join ('|',  map {"$_->[0]:$_->[1]"} @{$_[0]});
}
#_________________________________________________________________________________________

#	expand_cdna_with_gaps

#		input:

#			$seq
#			$aligning character
#			@list of position/length pairs in string format 0:63|119:6|263:1

#		returns:

#			$expanded sequences
#_________________________________________________________________________________________
sub expand_cdna_with_gaps($$$)
{
	my ($seq, $char, $align) = @_;
	$char |= '-';
	$char = substr($char, 0,  1);
	$char = $char x 3;
	
	foreach (reverse @$align)
	{
		substr($seq, $_->[0] * 3, 0) = ($char x $_->[1]);
	}
	return $seq;

}

#_________________________________________________________________________________________

#	expand_prot_with_gaps

#		input:

#			$seq
#			$aligning character
#			@list of position/length pairs

#		returns:

#			$expanded sequences
#_________________________________________________________________________________________
sub expand_prot_with_gaps($$\@)
{
	my ($seq, $char, $align) = @_;
	$char |= '-';
	$char = substr($char, 0,  1);
	
=pod
		#debug		
		print STDERR "length = ", length $seq, "\n";
		print STDERR "$_->[0],$_->[1];" for (@$align);
		print STDERR "\n";
=cut
	foreach (reverse @$align)
	{
		substr($seq, $_->[0], 0) = ($char x $_->[1]);
	}
	return $seq;

}

sub expand_cdna_with_gapstr($$$)
{
	my ($seq, $char, $gapstr) = @_;
	return expand_cdna_with_gaps($seq, $char, @{gapstr_to_gaps($gapstr)})
}

sub expand_prot_with_gapstr($$$)
{
	my ($seq, $char, $gapstr) = @_;
	return expand_prot_with_gaps($seq, $char, @{gapstr_to_gaps($gapstr)})
}



#_________________________________________________________________________________________

#	map_pos_to_aligned

#		input:

#			$position
#			@list of alignment position/length pairs

#		returns:

#			adjusted position given the aligning marks
#_________________________________________________________________________________________
sub map_pos_to_aligned($\@)
{
	my ($pos, $align) = @_;
	my $new_pos = $pos;
	foreach (@$align)
	{
		$new_pos += $_->[1] if ($pos >= $_->[0]);
	}
	return $new_pos;

}

#_________________________________________________________________________________________

#	map_pos_from_aligned

#		input:

#			$position
#			@list of alignment position/length pairs

#		returns:

#			adjusted position given the aligning marks
#_________________________________________________________________________________________
sub max($$)
{
	return $_[0] > $_[1] ? $_[0] : $_[1];
}
sub map_pos_from_aligned($\@)
{
	my ($pos, $gaps) = @_;

	# in gap region
	#if (scalar @_ && $gaps->[0][0] == $pos)
	#{
	#	return $gaps->[0][0] - 1;
	#}
	foreach (@$gaps)
	{
		# before next gap
		if ($_->[0] >= $pos)
		{
			if ($_->[0] == $pos)
			{
				return max($pos - 1, 0);
			}
			last;
		}
		
		$pos -= $_->[1];

		# in gap region
		if ($_->[0] > $pos)
		{
			$pos = $_->[0] - 1;
			last;
		}
	}
	return max($pos, 0);
}

#_________________________________________________________________________________________

#	map_pos_via_aligned


#		map from one alignment to another


#		input:

#			$position
#			@list of alignment position/length pairs for my sequence
#			@list of alignment position/length pairs for other sequence

#		returns:

#			adjusted position given the aligning marks
#_________________________________________________________________________________________
sub map_pos_via_aligned($\@\@)
{
	my ($pos,  $gaps1,  $gaps2) = @_;
	return map_pos_from_aligned(map_pos_to_aligned($pos, @$gaps1),  @$gaps2);
}













#_________________________________________________________________________________________

#	aligned_seq_to_gaps

#		input:

#			$ aligned protein sequence with gaps

#		returns:

#			positions of aligning marks
#_________________________________________________________________________________________
sub aligned_seq_to_gaps
{
	my ($seq, $gap_char) = @_;
	croak ("No sequence supplied") unless defined $seq;
	
	$gap_char = '-' unless defined $gap_char;
	
	my @gaps;
	my $cummulative = 0;
	# save the number and position of the aligning characters for each sequence
	while ($seq =~ /($gap_char+)/g)
	{
		# save position and length of aligning markers
		my $len = length($1);
		push(@gaps, [pos($seq) - $len - $cummulative, $len]);
		$cummulative += $len;
	}
	return [@gaps];
}



#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
#
#
#	Test
#

sub test
{
	my $PRINT_LEN = 25;
	my @align1 = ([11,3]);
	my @align2 = ([0,4],[11,3],[14,3]);
	my $seq1 = 'A' x 21;
	my $seq2 = 'B' x 14;
	
	expand_prot_with_gapstr($seq1,  '-',  '');
	expand_prot_with_gapstr($seq2,  '-',  '');
	
	
	#
	#	Sanity checks
	#
	for (0..25)
	{
		die unless map_pos_from_aligned(map_pos_to_aligned($_,  @align1), @align1) == $_;
		die unless map_pos_from_aligned(map_pos_to_aligned($_,  @align2), @align2) == $_;
	}
	die unless gaps_to_gapstr(@{aligned_seq_to_gaps(expand_prot_with_gaps($seq1, '-', @align1))}) eq 
				gaps_to_gapstr(@align1);
	die unless gaps_to_gapstr(@{aligned_seq_to_gaps(expand_prot_with_gaps($seq2, '-', @align2))}) eq 
				gaps_to_gapstr(@align2);
	
	
	#
	#	print conversion maps
	#
	print "\nseq1\t-> seq2\n----\t   ----\n";
	for (0..25)
	{
		print "$_\t-> ",
					map_pos_via_aligned($_, @align1, @align2), "\n";
	}
	
	
	my $expanded1 = expand_prot_with_gaps($seq1, '-', @align1);
	my $expanded2 = expand_prot_with_gaps($seq2, '-', @align2);
	# convert to numbers
	sub number_nongap_positions($)
	{
		my ($seq) = @_;
		my $counter = 0;
		my $numstr = '';
		for (split //, $seq)
		{
			if ($_ eq '-')
			{
				$numstr .= ' ';
			}
			else
			{
				$numstr .= ($counter % 10);
				++$counter;
			}
		}
		return $numstr;
	}
	
	print "\nSequences:\n";
	print substr(number_nongap_positions($expanded1), 0, $PRINT_LEN), "\n";
	print substr($expanded1, 0, $PRINT_LEN), "\n";
	print $_ % 10 for (0..$PRINT_LEN -1);
	print "\n";
	print substr($expanded2, 0, $PRINT_LEN), "\n";
	print substr(number_nongap_positions($expanded2), 0, $PRINT_LEN), "\n";
	print "\n";
}

#test;
1;

